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Abstract. This note addresses possible applications of the Tikhonov 
regularization to image reconstruction of gravitational lens systems. Sev- 
eral modifications of the regularization algorithm are discussed. Our illus- 
trative example is the close quadruple gravitational lens QSO 2237+0305 
(Einstein Cross). The restored image of the lens is decomposed into two 
parts - the quasar components and the background galaxy. 



1. Introduction 



Over the last decade gravitational lensing has become a powerful tool to probe 
the Universe. The number of lenses discovered so far is just below one hundred 
and is growing quickly. Gravitational lensing comes in a variety of shapes and 
sizes. One of the most interesting examples of this phenomenon is systems of 
close multiple quasar images. For a typical lens galaxy the image separation is 
of order of an arc second which happens to be the resolution of some of the best 
ground-based telscopes. Hence, what the observer will see is the source copies 
blurred and overlapped and possibly superimposed on the lens galaxy. To use 
the "nature telescope" one needs to restore the light distribution preferably 
splitting it into parts belonging to the "background" and to the quasar images. 
The problem addressed in this note is the problem of the correct reconstruction 
of complex gravitational lens images. 
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2. Algorithmic aspects of the image reconstruction 

There is a number of image reconstruction algorithms available for astronomers 
at present. Let us briefly review some of them. For example, in the CLEAN 
method one proceeds by subtracting the quasar components from the image 
in order to obtain underlying galaxy. The process is repeated iteratively until 
the termination criterion is satisfied, which usually requires visual inspection 
of the gray-scale distribution of residuals. Another wide-spread method is the 
Maximum Entropy deconvolution. Here the solution not only fits the noisy data 
but also maximizes so-called 'entropy'. Essentially, the entropy term is a measure 
of the number of different features in the image. In some cases this method 
can lead to a significant improvement in quality, however the very definition 
of the entropy forces the resulting image to be smooth, which sometimes can 
be considered as a drawback. In both the CLEAN and the MEM algorithms 
some kind of a 'default' image has to be introduced. The approach developed by 
Magain, Courbin and Sohy (MCS) requires only the PSF of the final (restored) 
image to be set by hand. The MCS method proved to work well for the data on 
which the previous algorithms had failed to give the optimal results. However, 
problems might appear when the MCS algorithm is faced with the image with 
irregular or varying over the field point spread function (PSF). 

2.1. Deconvolution basics 

The mathematical model of the image corruption due to the finite resolving 
power of the telescope and atmospheric perturbation is a Fredholm type integral 
equation of the first kind with a space-invariant kernel, or a convolution integral 
equation: 



Here u(x, y) represents the observed light distribution, t(x, y) is the point spread 
function determined from observations of the individual star at the periphery of 
the frame, z(x,y) represents the object, or sought solution, B is the frame area, 
B = [0, L] x [0, L\. The convolution operator A is the linear operator which 
acts from some Hilbert space Z to the Hilbert space of second-power-integrable 
functions L<i- The goal is to find an approximate solution having at our disposal 
the noisy data ug and the estimate of the noise level 5: \\u — us\\l 2 < This 
inverse problem belongs to the general class of ill-posed problems in the sense 
of Hadamar. 

There is a well-developed mathematical approach to solve such problems 
based on the idea of the regularizing algorithm. Since the publication of the 
fundamental work by A. Tikhonov (1963) this method has been extensively 
studied and widely adopted in many fields of science. 

2.2. A regularizing algorithm 

To obtain stable and physically valid results in image reconstruction the reg- 
ularization ought to be based on a priori knowledge of the properties of the 
admissible solution. 
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The regularization implies construction of the algorithm that controls the 
trade-off between the assumptions about smoothness and the structure of the 
sought solution and its consistency with the data. The key concept of the algo- 
rithm is a smoothing function: 

M a [z] = \\A[z] - us\\l 2 + a ■ n[z] (2) 

Here the first term represents the squared discrepancy, a is the regularization 
parameter, Q[z] is a stabilizer function. Let z a be the extremum of the function 
M a [z] on Z, i.e. z a is the solution of the minimization problem for M a [z] on 
the chosen set of functions (possibly with some constraints). 

The choice of the regularization parameter a is crucial for solving ill-posed 
problems. Generally, it should depend on the input data, the errors, and the 
method of approximation of the initial problem. One of the way to co-ordinate 
the regularization parameter with the error of the input information is the dis- 
crepancy principle: 

\\A[z a ]-u 5 \\ L2 ~6. (3) 

Providing the regularization parameter a is chosen according to this rule, z a 
can be considered as an approximate solution which tends to the exact solution 
in the context of the norm of the chosen set of functions as the error level of 
input data tends to zero. 

A prior knowledge about the smoothness of the unknown solution is embed- 
ded in the regularizing algorithm through the appropriate choice of the stabilizer 
function. In most cases it is the squared norm of the solution on some set of 
functions: Q[z] = \\z\\ 2 z . The choice of the stabilizer affects the order of con- 
vergence of approximate solutions. Thus, if it is assumed that the unknown 
solution belongs to the class of second-power-integrable functions, Z = L2, the 
stabilizer can be chosen as follows: 

to[z] = |Ml! 2 =JI z 2 dxdy, (4) 

B 

Regularization with this stabilizer function guarantees that approximate solu- 
tions converge to the exact solution in the context of the norm L2, i.e. converges 
in mean (zero-order convergence): 

Ik" — z\\l 2 ^0 as 5^0. 

If a priori information about sought solution allows to assume higher smooth- 
ness of z and choose Z = W22, where W22 is a set of L2-functions which have 
generalized derivatives of the second order those are second-power-integrable, 
the stabilizer can be written in the following form: 

o^,i*^//{^(£) 2 +2 (||) 2 + (0) 2 }^, ,5, 

When selecting a in accordance with the discrepancy principle, approximate 
solutions z a tend to the exact solution of the problem as 5 tends to zero in the 
context of the W22 norm: 

\\z a — z \\w 2 2 ~ * as $ ~~ * 0- 
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According to Sobolev's embedding theorem, 1^22 is embedded in C[B] 
the set of continuous functions on B. Thus, the convergence in the context 
of the W22-norm means the convergence in the context of the norm of C[B], i.e. 
regularized solutions converge to the exact solution uniformly: 

max \z a (x, y) — z(x, y)\ — > as 5^0 

(x,y)eB 

The degree of smoothness of the sought solution is not the only a priori 
information that can be considered and included into the regularizing algorithm. 
Various assumptions about the structure of the object under study can also be 
taken into account. Images of close quadruple gravitational lens systems consist 
of multiple overlapped quasar images superimposed on a background galaxy. 
So the image can be decomposed into two constituent parts - the sum of K 
5-functions and smooth background (galaxy): 

K 

z (x, V) = ak ^ x ~ h k,y - Cfe) + g(x, y), (6) 
fc=i 

where K is the number of point sources with coordinates (bk, C&) and intensities 
cifc in the frame; g(x, y) is the solution's component corresponding to a galaxy; 
5 represents Dirac function. 

Rapid intensity variations in the observed data caused by the bright nucleus 
of the galaxy can be processed by incorporation of an additional 5-function for 
the nucleus or by selecting the function describing the light distribution of the 
background galaxy from the appropriate set of functions, i.e. bounded total 
variation (TV) functions. 

An approach of piecewise uniform regularization based on TV class of func- 
tions was suggested by A. Leonov (1999). Let us consider an arbitrary grid 
SniN 2 introduced on B and define the total variation for a function z on B as 
follows: 

JVi-l JV 2 -1 
V(Z,B)= SUP ( £ |z m+ l,l - Z m:1 \ + \zi,n+l ~ Zl,n\ + 

Sn 1 n 2 171=1 n=1 

JVi-1 N 2 -l 

+ E E I Zm+l,n+l ^m+l,n ^m,n+l Z m ,n | j ^''S'A r iA r 2 ) 
m=l n=l 

The function for which the total variation is a finite quantity is called bounded 
total variation function. It is continious nearly everywhere with the exception, 
possibly, of the points of discontinuity positioned on the countable set of grid- 
lines. The regularization algorithm with the proper choice of the regularization 
parameter and the stabilizer function 

tt[z] = \\ z \\u[B] = \z(0,0)\+V(z,B) (7) 

provides piecewise uniform convergence of approximate solutions. 

Additionally, one can penalize the unknown solution for being drastically 
different from the certain analytical model and construct the stabilizer in the 
following form: 

M[z] = \\g - gmodelWc ( 8 ) 
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In this work we assume that the light distribution in the central region of the 
galaxy is well-modeled by generalized de Vaucouleurs profile (Sersic's model): 

9model{r) = Ioexp~ bn (~ )n (9) 
where b n = 2n - 0.324 for 1 < n < 4. 



3. Numerical results 



3.1. Observations 



The observations of QSO 2237+0305 were carried out during August, September 
and October of 2002 on the 1.5-m telescope AZT-22 at the Maidanak Obser- 
vatory in Uzbekistan. The CCD camera with the gain of g = 1.2e~ ADU^ 1 , 
the readout noise of RON = 10e~ and the pixel scale of Q.12"pixel~ l was used. 
In this work one of R-band frames obtained on the 29th of August 2002 was 
taken. The best quality of the image corresponds to the point source with the 
FWHM=0.75". Preprocessing of the data including bias-level subtraction, flat- 
field division, sky subtraction and cosmic ray removal was done with the stan- 
dard routines in Munich Image Data Analysis System (MIDAS) environment. 
The set of data taken over the night was averaged and the constant background 
over the frame was subtracted to prepare image of Einstein Cross for the further 
treatment. Then the subframe of 64 by 64 pixels with the image of Einstein 
Cross centered on nucleus of the galaxy 2237+0305 was extracted. The total 
noise over the frame was calculated as follows: 



\ ~f Ng N K ' 

where N=4 is the number of the frame of one observation night; cross^ - intensity 
of ij-pixel. For the subframe the total noise is about 400 ADU. 

3.2. PSF processing 

In the Einstein Cross data, the quasar images lie so close to each other that 
the PSF from each image overlap with the PSFs of its neighbours. But for 
normal atmospheric conditions it is vital to have a clean isolated PSF for the 
data frame that we want to process. The only way then is to derive the PSF 
from the bright star labelled a in the frame. To create the PSF the star was 
extracted from the frame. It has to be noted that the shape of the PSF can 
be very complicated. This includes the sum of the imaging properties of the 
telescope and the effects of atmospheric turbulence on the signal. There is a 
number of methods to construct a PSF. In this work both numerical PSF and 
Gaussian PSF were tested. The use of a numerical PSF gave the best fit to 
the quasar components. The wings of the PSF light distribution were slightly 
smoothed by the median filter with the window size compared to the FWHM. 
To find the parameters of Gaussian PSF it is necessary to solve the inverse 
problem. However, the approximation of the PSF by the Gauss profile is not 
a good representation of the real light distribution. If k is the ratio of total 
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Figure 1. PSF fitting: a) numerical PSF (star from frame); b) resid- 
uals (GaussianPSF — numerical PSF) / V GaussianPSF; c) residuals 
(filteredPSF - numericalPSF) /^filteredPSF 



intensity of the model and the total intensity of the observed star, the value of 
summarised intensity losses are derived as follows: 



y = (1 - k)100% (11) 

In the case of the Gaussian PSF the losses amount to 4.8% and in the case of 
numerical PSF slightly smoothed by filter the losses do not exceed 3.9%. 



3.3. Image reconstruction for QSO 2237+0305 

Working with CCD-frames one has a pixel grid naturally introduced in the image 
frame and deals with discrete functions in (1). So it is necessary to consider the 
finite difference approximation of the smoothing function: 



N -i \ N ( / 4 



2 



M a [a k , gij ] = ~r { J2 { h-hg-j a ^i-b k ,j-c k + 9ij )\ - u pq \ +uQ.[ gij ] 



p,q=l a PQ U,j=l k \fe=l 



(12) 

Here the variables are the coordinates and the intensities of the quasar compo- 
nents and the pixel values of the background galaxy. 



Image reconstruction on L 2 The regularization on L 2 set of functions repre- 
sents the simplest sort of Tikhonov's method. Written for the pixel grid the 
stabilizer: 

N 

m^] = E 9% (13) 

Here g = g — g S ersic- The parameters of the Sersic model were obtained at 
preliminary stage using the least-squares method. For the minimization of the 
smoothing function (12) the conjugate- gradient method was used. Figure 2 
shows the results of the image reconstruction. The map of residuals was calcu- 
lated as follows: (modelij — crosstj) / aij. Figure 3 shows the astrometry results. 
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Figure 2. Reconstruction on L<i'- a) Einstein Cross data; b) solution 
on Z/2 (logarithmic scale); c) solution convolved with PSF; d) residuals 
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Figure 3. Astrometry for Z/2-reconstruction 
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Figure 4. Reconstruction on W22: a) the Einstein Cross data; b) 
solution on W22 (logarithmic scale); c) solution convolved with PSF; 
d) residuals 
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Figure 5. Astrometry for Vl^-reconstruction 



Image reconstruction on W22 The discrete expression for the stabilizer is: 

N N-l N 

^[9ij] = E 9ij + E E - 2 9i,j + k-i,j) 2 + 

ij=l i=2 j=l 

N-l N-l N N-l 

E E 2(^+1^+1 - gi+ij - g iij+ i + gtj) 2 + E E (di,j+i ~ 2 9i,j + 9i,j-i) 2 

i=l j=l i=l j=2 

Figure 4 shows the results of the image reconstruction with the regularizing 
algorithm based on the W22 set of functions. 

Image reconstruction on v[B] The regularization technique with the use of 
bounded total variation functions is applicable if there are reasons to assume 
that the unknown solution has discontinuities or rapid variations. The discrete 
expression for stabilizer reads: 

N-l 

^l z ] = fe(9i+l,j+l ~ 9i+i,j - 9i,j+l + 9i,j) 
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a] 




Figure 6. Reconstruction on the v[B\: a) Einstein Cross data; b) 
solution on v[B\; c) solution convolved with PSF; d) residuals 



where f e (t) = yt 2 + yj^ij- The results of the image restoration are presented 
on Figure 6. 

4. Conclusions 

The application of the Tikhonov regularization to the image reconstruction of 
gravitational lens systems seems to be an effective way of tackling some of the 
common problems shared by conventional deconvolution algorithms. We have 
considered several modifications of the regularization method. The illustrative 
example of the Einstein Cross lens system was used to compare different versions 
of the algorithm. Although the results look quite encouraging the further work 
needs to be done to select the appropriate stabilizer function and to produce the 
pipe-line version of the software. 
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